#' plot.bapvarIRF
#'
#' Plot IRFs from Posterior BaP-VAR Draws
#'
#' @param x An object of class \code{bapvarIRF}
#'   (see \code{\link{bapvar_dynamics}})
#' @param var_names A character vector giving the endogenous variable names
#'   (optional; if missing, plots labeled as "Equation 1", "Equation 2", ...)
#' @param line_color Colors to draw the IRF estimate lines
#'   (optional; if missing, the lines are drawn using colors generated by
#'   \code{\link[bggum]{okabe_ito}()})
#' @param ci_color Colors to shade the credible intervals (optional;
#'   if missing, it is shaded using colors generated by
#'   \code{\link[bggum]{okabe_ito}()}, with alpha value 37 (of 255)).
#' @param ci A numeric vector of length one giving the credible interval
#'   desired; default is 0.68
#' @param est A character vector of length one; should the estimate use the
#'   mean or the median? (default is "median")
#' @param na.rm A logical vector of length one; should NA values be removed?
#'   (default is TRUE)
#' @param reorder A logical vector of length one; should any ordering done
#'   by \code{bapvar_irf} by undone? (default is TRUE)
#' @param ... Arguments passed to or from other methods (currently unused)
#'
#' @return NULL (invisibly)
#' @export
plot.bapvarIRF <- function(x, var_names, line_color, ci_color,
                           ci = 0.68, est = "median",
                           na.rm = TRUE, reorder = TRUE, ...) {
    if ( ! "bapvarIRF" %in% class(x) ) {
        stop("Provide a bapvarIRF object. See help(\"bapvar_dynamics\".")
    }
    # Get ordering
    ordering <- attr(x, "ordering")
    # Summarize IRF draws
    qs <- c( (1 - ci) / 2, 0.5, (1 + ci) / 2 )
    x <- apply(x, 2:3, quantile, probs = qs, na.rm = na.rm)
    if ( est == "mean " ) {
        x[2, , ] <- apply(x, 2:3, mean, na.rm = na.rm)
    }
    # Record dimensions
    dims <- dim(x)
    m <- sqrt(dims[2])
    h <- dims[3]
    # Deal with variable names and ordering
    if ( missing(var_names) ) {
        var_names <- paste("Equation", 1:m)
    }
    if ( missing(line_color) ) {
        line_color <- bggum::okabe_ito(m)
    }
    if ( missing(ci_color) ) {
        ci_color <- paste0(bggum::okabe_ito(m), "37")
    }
    if ( reorder ) {
        mat_ordering <- rep(ordering, m) + rep(ordering-1, each = m) * m
        mat_ordering_original <- 1:(m^2)
        new_ordering <- match(mat_ordering_original, mat_ordering)
        x <- x[ , new_ordering, ]
    }
    # Work out margins + "layout"
    opar <- par(mfcol = c(m, m), mar = c(2, 2, 1, 1) + 0.1, oma = c(1, 5, 6, 0))
    on.exit(par(opar))
    # Alignments for axis labels
    # XXX: This really only works for m such that 2 <= m <= 10;
    # I can't envision that ever being an issue, but worth noting
    k <- floor(m/2)
    adjustments <- c(-0.05 / (1:k), "if"(m %% 2 == 0, NULL, 0), 0.05 / (k:1))
    alignments <- (c(0, (1:(m-1)/m)) + (1:m)/m) / 2 + adjustments
    # Make plots
    for ( j in 1:m ) {
        for ( i in 1:m ) {
            lows <- x[1, (j-1)*m + i, ]
            medians <- x[2, (j-1)*m + i, ]
            highs <- x[3, (j-1)*m + i, ]
            lims <- range(c(lows, medians, highs))
            plot(medians, type = "n", ylim = lims, xlab = "", ylab = "",
                 xaxt = "n", yaxt = "n", main = "")
            polygon_x <- c(1:h, h:1)
            polgyon_y <- c(lows, rev(highs))
            polygon(polygon_x, polgyon_y, col = ci_color[i], border = NA)
            abline(h = 0, lty = 2, col = "black")
            lines(1:h, medians, type = "o", pch = 19, col = line_color[i])
            axis(1, at = 1:h, labels = 0:(h-1), tick = FALSE, line = -0.75,
                 cex.axis = 1.5)
            axis(2, tick = FALSE, line = -0.75, cex.axis = 1.5)
            if ( j == 1 ) {
                mtext(var_names[i], side = 2, line = 0.5, cex = 1.5,
                      outer = TRUE, adj = alignments[m + 1 - i])
                if ( i == 2 ) {
                    mtext("Shock to", side = 2, line = 3, cex = 1.5,
                          outer = TRUE)
                }
            }
            if ( i == 1 ) {
                mtext(var_names[j], line = 0.5, cex = 1.5, outer = TRUE,
                      adj = alignments[j])
                if ( j == 2 ) {
                    mtext("Response in", line = 3, cex = 1.5, outer = TRUE)
                }
            }
        }
    }
    invisible(NULL)
}

#' plot.bapvarFEVD
#'
#' Plot FEVDs from Posterior BaP-VAR Draws
#'
#' @param x An object of class \code{bapvarFEVD}
#'   (see \code{\link{bapvar_dynamics}})
#' @param var_names A character vector giving the endogenous variable names
#' @param color_palette A vector of colors to plot the bars;
#'   if missing, the bars are filled using colors generated by
#'   \code{\link[bggum]{okabe_ito}()}
#' @param na.rm A logical vector of length one; should NA values be removed?
#'   (default is TRUE)
#' @param reorder A logical vector of length one; should any ordering done
#'   by \code{bapvar_irf} by undone? (default is TRUE)
#' @param ... Arguments passed to or from other methods (currently unused)
#'
#' @return NULL (invisibly)
#' @export
plot.bapvarFEVD <- function(x, var_names, color_palette,
                            na.rm = TRUE, reorder = TRUE, ...) {
    if ( ! "bapvarFEVD" %in% class(x) ) {
        stop("Provide a bapvarFEVD object. See help(\"bapvar_dynamics\".")
    }
    # Record dimensions
    dims <- dim(x)
    m <- sqrt(dims[2])
    h <- dims[3]
    ordering <- attr(x, "ordering")
    # Summarize FEVDs
    x <- apply(x, 2:3, mean, na.rm = na.rm)
    # Reorder (or not)
    if ( reorder ) {
        mat_ordering <- rep(ordering, m) + rep(ordering-1, each = m) * m
        mat_ordering_original <- 1:(m^2)
        new_ordering <- match(mat_ordering_original, mat_ordering)
        x <- x[new_ordering, ]
    }
    # Set up titles, colors, and "layout"
    if ( missing(var_names) ) {
        var_names <- paste("Equation", 1:m)
    }
    if ( missing(color_palette) ) {
        color_palette <- bggum::okabe_ito(m)
    }
    if ( length(color_palette != m) ) {
        color_palette <- rep(color_palette, length.out = m)
    }
    plot_titles <- paste("Response in", var_names)
    legend_text <- paste("Shock to", var_names, "  ")
    opar <- par(mfcol = c(1, m), mar = c(3, 3, 3, 1) + 0.1, oma = c(3, 3, 0, 0))
    on.exit(par(opar))
    middle_eqn <- ceiling(m / 2)
    xjust <- "if"(middle_eqn %% 2 == 0, 0.5, 0)
    adjust <- "if"(middle_eqn %% 2 == 0, 0, -h/4)
    # Plot FEVDs
    for ( j in 1:m ) {
        these_fevds <- x[1:m + (j-1)*m, ]
        mp <- barplot(these_fevds, col = color_palette, xlab = "", ylab = "",
                      main = plot_titles[j], xaxt = "n", yaxt = "n",
                      ylim = c(0, 1), cex.main = 1.5)
        xat <- mp
        xlab <- 0:(h-1)
        yat <- seq(from = 0, to = 1, by = 0.25)
        ylab <- c("0", "0.25", "0.5", "0.75", "1")
        axis(side = 1, at = xat, labels = xlab, tick = FALSE, line = -0.75,
             cex.axis = 1.5)
        axis(side = 2, at = yat, labels = ylab, tick = FALSE, line = -0.75,
             cex.axis = 1.5)
        if ( j == 1 ) {
            mtext("FEVD", side = 2, line = 0.5, cex = 1.5, outer = TRUE)
        }
        if ( j == middle_eqn ) {
            X <- mean(mp) + adjust
            legend(x = X, y = -0.1, legend = legend_text, fill = color_palette,
                   bty = "n", horiz = TRUE, xpd = NA, cex = 2, xjust = xjust)
        }
    }
    invisible(NULL)
}
